* Please ensure the working directory is the same as
* that for selectionModel.do (this file)

* Sample selection model for willingness to pay


*********************************
* Import .csv data
*********************************

* If generating data from R file, note that you need to remove "NA" cells or you will receive a "type mismatch" error on numeric values.

import delimited "first_stage.csv", clear


*********************************
* Set up data
*********************************

gen lackAccess = 1 if use_grid == "No"
replace lackAccess = 0 if use_grid == "Yes"

gen willPay = willingness_to_pay if lackAccess == 1

gen bc = 1 if unit_caste != "General"
replace bc = 0 if unit_caste == "General"

gen educ10 = 1 if unit_educ == "Up to 10th Standard"
replace educ10 = 1 if unit_educ == "12th Standard or Diploma"
replace educ10 = 1 if unit_educ == "Graduate and Above"
replace educ10 = 0 if educ10 == . & unit_educ != ""

gen lwillPay = ln(willPay + 1)

gen lham_hh = ln(ham_hh)

gen lcens_distance_town = ln(cens_distance_town + 1)

gen lunit_expend = ln(unit_expend)

encode state, gen(nstate)

tab state, gen(state)

sort village_code
by village_code: egen m_pr_kerosene_mkt = mean(pr_kerosene_mkt)
by village_code: egen m_pr_kerosene_pds = mean(pr_kerosene_pds)

gen kerosene_expend3 = (kerosene_mkt * pr_kerosene_mkt) + (kerosene_pds * pr_kerosene_pds)
replace kerosene_expend3 = (kerosene_mkt * m_pr_kerosene_mkt) + (kerosene_pds * m_pr_kerosene_pds) if kerosene_expend3 == . & kerosene_mkt < .
replace kerosene_expend3 = (kerosene_mkt * m_pr_kerosene_mkt) + (kerosene_pds * m_pr_kerosene_pds) if kerosene_expend3 == . & kerosene_pds < .
replace kerosene_expend3 = 0 if kerosene_expend3 == .
gen ln_kerosene_expend3 = ln(kerosene_expend3 + 1)


***************************************
* OLS Regresssion Analyses
***************************************

* Analysis of overall hours

eststo clear

eststo: reg lwillPay ev_hrs_d lham_hh lcens_distance_town lunit_expend ///
  educ10 bc ln_kerosene_expend3 i.nstate, vce(cluster village_code)
estadd local state "Yes" , replace

margins, at(ev_hrs_d=(1 3 5 7 9 11 13 15 17 19 21 23))
marginsplot, ytitle(ln(Willingness to Pay)) ylabel(0(2)8) xtitle(Available Hours) scheme(s1mono) name(plot1)

* Analysis of quality index

eststo: reg lwillPay ev_qual_d lham_hh lcens_distance_town lunit_expend ///
  educ10 bc ln_kerosene_expend3 i.nstate, vce(cluster village_code)
estadd local state "Yes" , replace

margins, at(ev_qual_d=(0 0.2 0.4 0.6 0.8 1))
marginsplot, ytitle(ln(Willingness to Pay)) ylabel(0(2)8) xtitle(Quality Index) scheme(s1mono) name(plot2)

* Analysis of nighttime hours

eststo: reg lwillPay ev_hrs_n lham_hh lcens_distance_town lunit_expend ///
  educ10 bc ln_kerosene_expend3 i.nstate, vce(cluster village_code)
estadd local state "Yes" , replace

margins, at(ev_hrs_n=(0 1 2 3 4 5 6))
marginsplot, ytitle(ln(Willingness to Pay)) ylabel(0(2)8) xtitle(Nighttime Hours) scheme(s1mono) name(plot3)


gr combine plot1 plot2 plot3, scheme(s1mono)
graph drop _all


* Create table

esttab using "Table_4.tex", ///
   replace b(%9.3f)  ///
   stats(state vce r2 N, fmt(%~12s %~12s %9.3f %9.3g) ///
   label("Fixed effects: state" "VCE" "R-squared" "Observations")) ///
   varlabels(ev_hrs_d "Available Hours" ev_qual_d ///
     "Quality Index" ev_hrs_n "Nighttime Hours" lham_hh ///
     "ln(Habitation Households)" lcens_distance_town "ln(Distance to Town)" ///
     lunit_expend "ln(Household Expenditure)" educ10 "10th Year or Greater Education" ///
     bc "Backwards Caste" ln_kerosene_expend3 "ln(Kerosene Expenditure)" _cons Constant) ///
   booktabs eqlabels(none) ///
   drop(1.nstate 2.nstate 3.nstate 4.nstate 5.nstate 6.nstate) ///
   order(ev_hrs_d ev_qual_d ev_hrs_n) ///
   se label star(* 0.10 ** 0.05 *** 0.01) nomtitles
 

***************************************
* Selection Analyses
***************************************

eststo clear

* Selection model for overall hours

eststo: heckman lwillPay ev_hrs_d lcens_distance_town lunit_expend educ10 bc ln_kerosene_expend3 ///
   i.nstate, select(lackAccess = lham_hh ev_hrs_d lcens_distance_town ///
   lunit_expend educ10 bc ln_kerosene_expend3 i.nstate) vce(cluster village_code)
estadd local state "Yes" , replace

margins, at(ev_hrs_d=(1 3 5 7 9 11 13 15 17 19 21 23))
* Note the margins are off because State does not yet include 1st stage calculations in 
* their implementation of margins. These are manually deflated in charts.
predict hdoos in 6456/6467, ycond


* Selection model for quality index

eststo: heckman lwillPay ev_qual_d lcens_distance_town lunit_expend educ10 bc ln_kerosene_expend3 ///
   i.nstate, select(lackAccess = lham_hh ev_qual_d lcens_distance_town ///
   lunit_expend educ10 bc ln_kerosene_expend3 i.nstate) vce(cluster village_code)
estadd local state "Yes" , replace

margins, at(ev_qual_d=(0 0.2 0.4 0.6 0.8 1))
predict qloos in 6456/6467, ycond
 
* Selection model for nighttime hours

eststo: heckman lwillPay ev_hrs_n lcens_distance_town lunit_expend educ10 bc ln_kerosene_expend3 ///
   i.nstate, select(lackAccess = lham_hh ev_hrs_n lcens_distance_town ///
   lunit_expend educ10 bc ln_kerosene_expend3 i.nstate) vce(cluster village_code)
estadd local state "Yes" , replace

margins, at(ev_hrs_n=(0 1 2 3 4 5 6))
predict hnoos in 6456/6467, ycond

gr combine plot1 plot2 plot3, scheme(s1mono)
graph drop _all


*Create Table

esttab using "Table_5.tex", ///
   replace b(%9.3f)  ///
   stats(state vce N, fmt(%~12s %~12s %9.3g) ///
   label("Fixed effects: state" "VCE" "R-squared" "Observations")) ///
   varlabels(ev_hrs_d "Available Hours" ev_qual_d ///
     "Quality Index" ev_hrs_n "Nighttime Hours" lham_hh ///
     "ln(Habitation Households)" lcens_distance_town "ln(Distance to Town)" ///
     lunit_expend "ln(Household Expenditure)" educ10 "10th Year or Greater Education" ///
     bc "Backwards Caste" ln_kerosene_expend3 "ln(Kerosene Expenditure)" _cons Constant) ///
   booktabs eqlabels(none) ///
   drop(1.nstate 2.nstate 3.nstate 4.nstate 5.nstate 6.nstate) ///
   order(ev_hrs_d ev_qual_d ev_hrs_n) ///
   se label star(* 0.10 ** 0.05 *** 0.01) nomtitles

eststo clear

************************************************
**Check relationship between habitation size and proportion of households running businesses
************************************************

cor ham_business lham_hh

eststo clear

* Selection model for overall hours

eststo: heckman lwillPay ev_hrs_d lcens_distance_town lunit_expend educ10 bc ln_kerosene_expend3 ham_business ///
   i.nstate, select(lackAccess = lham_hh ev_hrs_d lcens_distance_town ///
   lunit_expend educ10 bc ln_kerosene_expend3 ham_business i.nstate) vce(cluster village_code)
estadd local state "Yes" , replace


* Selection model for quality index

eststo: heckman lwillPay ev_qual_d lcens_distance_town lunit_expend educ10 bc ln_kerosene_expend3 ham_business ///
   i.nstate, select(lackAccess = lham_hh ev_qual_d lcens_distance_town ///
   lunit_expend educ10 bc ln_kerosene_expend3 ham_business i.nstate) vce(cluster village_code)
estadd local state "Yes" , replace

* Selection model for nighttime hours

eststo: heckman lwillPay ev_hrs_n lcens_distance_town lunit_expend educ10 bc ln_kerosene_expend3 ham_business ///
   i.nstate, select(lackAccess = lham_hh ev_hrs_n lcens_distance_town ///
   lunit_expend educ10 bc ln_kerosene_expend3 ham_business i.nstate) vce(cluster village_code)
estadd local state "Yes" , replace

*Create Table

esttab using "Table_A1.tex", ///
   replace b(%9.3f)  ///
   stats(state vce N, fmt(%~12s %~12s %9.3g) ///
   label("Fixed effects: state" "VCE" "R-squared" "Observations")) ///
   varlabels(ev_hrs_d "Available Hours" ev_qual_d ///
     "Quality Index" ev_hrs_n "Nighttime Hours" lham_hh ///
     "ln(Habitation Households)" lcens_distance_town "ln(Distance to Town)" ///
     lunit_expend "ln(Household Expenditure)" educ10 "10th Year or Greater Education" ///
     bc "Backwards Caste" ln_kerosene_expend3 "ln(Kerosene Expenditure)" ham_business "Habitation Businesses" _cons Constant) ///
   booktabs eqlabels(none) ///
   drop(1.nstate 2.nstate 3.nstate 4.nstate 5.nstate 6.nstate) ///
   order(ev_hrs_d ev_qual_d ev_hrs_n) ///
   se label star(* 0.10 ** 0.05 *** 0.01) nomtitles

eststo clear
